Замечание: у меня начала заедать буква 'Б' на клавиатуре, поэтому я могла иногда ее пропустить и этого не заметить, ровно как и запятые. Я старалась все тщательно перечитать, но могла не все подправить. Так что заранее извиняюсь за такое написание!
Библиотеки¶
%load_ext autoreload
%autoreload 2
import os
import pandas as pd
import numpy as np
import zipfile
import shutil
from tqdm import tqdm
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("darkgrid", {"axes.facecolor": ".9"})
import warnings
warnings.simplefilter('ignore')
from task2.polysomy_search import KlinefelterSearcher
from task3.find_bam_assembly import BamAssemblyFinder
from task4.analyze_gene_cov import GeneCoverageAnalyzer
from task5.find_vcf_assembly import VCFAssemblyFinder
from task6.pair_compare_vcf import VCFComparator
from task7.analyze_genetic_variants import GeneticVariantAnalyzer
from task7.plot_3d_scatter import plot_3d_scatter
Решение¶
- Заходим в IGV и выбираем референсный геном hg19
- Загружаем BAM и BAI файлы (File $\rightarrow$ Load from URL...)
- В IGV в поле поиска вводим ген APOE - таким образом мы перемещаемся в нужную область генома
- Определяем аллели APOE. Для этого бежим в SNPедию и понимаем, что основные варианты гена APOE связаны с полиморфизмами rs429358 и rs7412
| rs429358 | rs7412 | Name |
|---|---|---|
| C | T | ε1 |
| T | T | ε2 |
| T | C | ε3 |
| C | C | ε4 |
Координаты на hg19 ищем в Varsome:
- Возвращаемся в IGV и сравниваем свои данные в этих позициях с референсным геномом:
- rs429358: на референсном геноме T. У нас гомозиготный вариант T/T
- rs7412: на референсном геноме C. У нас гомозиготный вариант C/C
- Теперь снова идем в SNPедию и оцениваем, какой риск Альцгеймера соответствует нашему генотипу (точнее генотипу нашего образца, а не нашему). Magnitude показывает степень риска, где 6 указывает на очень высокий риск, а 1 - на низкий риск развития АЛьцгеймера.
| Common name | Genoset | Magnitude | rs429358 | rs7412 | Comment |
|---|---|---|---|---|---|
| Apo-ε1/ε1 | gs267 | 6 | (C;C) | (T;T) | the rare missing allele |
| Apo-ε1/ε2 | gs271 | 2.5 | (C;T) | (T;T) | |
| Apo-ε1/ε3 | gs270 | 2.6 | (C;T) | (C;T) | ambiguous ε2/ε4 or ε1/ε3 |
| Apo-ε2/ε4 | gs270 | 2.6 | (C;T) | (C;T) | ambiguous ε2/ε4 or ε1/ε3 |
| Apo-ε1/ε4 | gs272 | 2.5 | (C;C) | (C;T) | |
| Apo-ε2/ε2 | gs268 | 4 | (T;T) | (T;T) | good; lowest risk |
| Apo-ε2/ε3 | gs269 | 2 | (T;T) | (C;T) | |
| Apo-ε3/ε3 | gs246 | 2 | (T;T) | (C;C) | the most common |
| Apo-ε3/ε4 | gs141 | 3 | (C;T) | (C;C) | |
| Apo-ε4/ε4 | gs216 | 6 | (C;C) | (C;C) | ~11x increased Alzheimer's risk |
Вывод
Для нашего образца характерен наиболее распространенный статус APOE: Apo-ε3/ε3, который указывает на низкий риск развития Альцгеймера (2 по 6-бальной шкале)
Решение¶
Идея¶
У мужчин с синдромом Клайнфельтера имеется как минимум одна лишняя X-хромосома (кариотип 47,XXY вместо нормального 46,XY). То есть у мужчин с данным синдромом 2 (или более) X-хромосомы.
У здорового мужчины количество прочтений, приходящихся на X-хромосому, должно быть примерно такое же, как на Y-хромосому. Это означает, что покрытие (coverage) обеих хромосом будет схожим, так как у мужчин по одной X и Y хромосоме. Однако при синдроме Клайнфельтера ожидается, что количество ридов на X-хромосоме будет примерно в 2 раза больше, чем на Y-хромосоме, так как у таких людей две X-хромосомы.
Ход решения¶
Извлечение количества ридов для каждой хромосомы. Для кадого из
.bamфайлов извлекаем количество ридов, пришедшихся на каждуюиз хромосом. Для этого будем использоватьsamtools. При этом,- $C_X$ - количество прочтений для X-хромосомы
- $C_Y$ - количество прочтений для Y-хромосомы
Нормализация. Нужно выполнить нормализацию на длину хромосомы. Когда мы говорим о количестве прочтений на участок, то есть сильная зависимость от длины этого участка. Чем он будет длиннее, тем больше ридов на него (скорее всего) придется. Поэтому корректнее будет сравнивать количество прочтений, нормализованное на длину каждой хромосомы. $$N_X = \frac{C_X}{L_X} \quad N_Y = \frac{C_Y}{L_Y},$$ где $L_X$ и $L_Y$ - длина $X$ и $Y$ хромосом соответственно
Ожидаем, что $\text{Ratio} =\frac{N_X}{N_Y}$ будет ранжироваться следующим образом: $$\text{Ratio}[\text{female}] < \text{Ratio}[\text{health\_male}] < \text{Ratio}[\text{disease\_male}]$$
- Визуализация. Построим столбчатую диаграмму для каждого образца, которая будет отражать $N_X$ и $N_Y$. Так как у различных образцов различная глубина секвенирования, то для более красивой и корректной визуализации будем нормализовать значения $N_X$ и $N_Y$ на $N_{total}$ (нормализованное суммарное число ридов)
# Указываем путь к папке с .bam и .bam.bai файлами
DATA_FOLDER_TASK2 = "task2/data/downsampled"
searcher = KlinefelterSearcher(data_folder=DATA_FOLDER_TASK2)
data_karyo = searcher.get_normalized_reads(save_data=True)
data_karyo.head(2)
Process .bam files: 100%|██████████| 26/26 [00:00<00:00, 141.78it/s]
| N_X | N_Y | N_total | N_X_norm | N_Y_norm | Ratio | |
|---|---|---|---|---|---|---|
| HG00105 | 0.000012 | 6.686477e-06 | 0.000021 | 0.556014 | 0.319458 | 1.740491 |
| HG00099 | 0.000024 | 2.357952e-07 | 0.000024 | 1.015929 | 0.009948 | 102.124853 |
searcher.plot_karyotypes(data_karyo, save_fig=True)
Выводы¶
- Образцы, $\text{Ratio}$ которых превышает десятки - принадлежат женскому полу (97, 99, 100, 102, 106, 1860). Мы видим какое-то число ридов, которые выровнялись на Y-хромосому, но это случайное событие. Такие риды пришли из других участков генома.
- Образцы, $\text{Ratio}$ которых меньше двойки - здоровые мужчины (101, 105, 107, 108, 109). Значения $\text{Ratio}$ у них не равняются ровно единице, потому что при нормировке мы не могли учесть все факторы, которые определяют число ридов, которое упадет на тот или иной участок (хромосому). Но с нашим априорным знанием о том, что большинство мужчин здоровы, и результатом в несколько похожих значений $\text{Ratio}$, мы можем заключить, что у здоровых мужчин $\text{Ratio}$ может быть как минимум в диапазоне [1.35, 1.75]
- Образец 1337 имеет $\text{Ratio} = 3.48$, что отличает его ото всех остальных образцов. Данное значение сильно выше, чем у мужчин (как мы решили, здоровых), но и ниже, чем у женщин. Так что у пациента 1337 кариотип 47,XXY и он обладатель синдрома Кляйнфейтера
Решение¶
В принципе в этом задании расписывать особо нечего - просто с помощью samtools достаем метаинформацию о .bam файле (header), где можно найти информацию о сборке файл. Дальше нужно "запарсить" некоторые теги, которые могут содержать информацию о сборке (что указано в документации samtools):
- @SQ - Reference sequence dictionary, в поле AS=Genome assembly identifier. Сразу скажу, что в этом файле я его не нашла, поэтому:
- @PG - Program. В общем, можно попытаться извернуться и просмотреть все команды, которыми получали
.bam-файл. Там должно быть название референсного генома, в котором может быть отсылка к сборке.
У нас данные в .zip-формате. Чтобы не плодить файлы и не занимать память, извлечем файлы во временную директорию и после удалим, оставив только изначально сжатый формат.
DATA_FOLDER_TASK3 = "task3/data"
task3_bampath = os.path.join(DATA_FOLDER_TASK3, "task3.bam.assembly.zip")
extracted_dir = os.path.join(DATA_FOLDER_TASK3, "tmp")
# 1. Распаковываем архив во временную папку
with zipfile.ZipFile(task3_bampath, 'r') as zip_ref:
zip_ref.extractall(extracted_dir)
# 2. Достаем из временной папки название .bam файла
for file_name in os.listdir(extracted_dir):
if file_name.endswith(".bam"):
bam_file = os.path.join(extracted_dir, file_name)
# Посмотрим на красивое название
bam_file
'task3/data/tmp/vi0006.genome.markdup.downsampled.bam'
# 3. Наконец-то найдем информацию о сборке
bam_assembly_finder = BamAssemblyFinder(bam_file_path=bam_file)
bam_assembly_finder.find_assembly_info()
Сборка (AS в теге @SQ) не найдена Найдена информация о программе выравнивания (@PG): @PG ID:bwa PN:bwa VN:0.7.17-r1188 CL:bwa mem -t 30 -K 10000000 -k 30 -T 30 -R @RG\tSM:vi0006\tID:genome\tPL:illumina /home/xsukhanova/ref/hg38/hg38.fna vi0006.genome.cut_1.fastq.gz vi0006.genome.cut_2.fastq.gz @PG ID:samtools PN:samtools PP:bwa VN:1.19.2 CL:samtools view -@ 30 -Sb - @PG ID:samtools.1 PN:samtools PP:samtools VN:1.19.2 CL:samtools sort -@ 30 -o vi0006.genome.sort.bam - @PG ID:MarkDuplicates VN:Version:4.2.0.0 CL:MarkDuplicates --INPUT vi0006.genome.sort.bam --OUTPUT vi0006.genome.hg38.markdup.bam --METRICS_FILE vi0006.genome.markdup.metrics.txt --MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP 50000 --MAX_FILE_HANDLES_FOR_READ_ENDS_MAP 8000 --SORTING_COLLECTION_SIZE_RATIO 0.25 --TAG_DUPLICATE_SET_MEMBERS false --REMOVE_SEQUENCING_DUPLICATES false --TAGGING_POLICY DontTag --CLEAR_DT true --DUPLEX_UMI false --ADD_PG_TAG_TO_READS true --REMOVE_DUPLICATES false --ASSUME_SORTED false --DUPLICATE_SCORING_STRATEGY SUM_OF_BASE_QUALITIES --PROGRAM_RECORD_ID MarkDuplicates --PROGRAM_GROUP_NAME MarkDuplicates --READ_NAME_REGEX <optimized capture of last three ':' separated fields as numeric values> --OPTICAL_DUPLICATE_PIXEL_DISTANCE 100 --MAX_OPTICAL_DUPLICATE_SET_SIZE 300000 --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_secrets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false PN:MarkDuplicates @PG ID:samtools.2 PN:samtools PP:samtools.1 VN:1.19.2 CL:samtools view -@ 45 -s 0.028 -bo vi0006.genome.hg38.markdup.downsampled.x1.bam vi0006.genome.hg38.markdup.bam @PG ID:samtools.3 PN:samtools PP:samtools.2 VN:1.20 CL:samtools view -h -b -o vi0006.genome.markdup.downsampled.bam vi0006.genome.markdup.downsampled.x1.bam chr20 @PG ID:samtools.4 PN:samtools PP:samtools.3 VN:1.20 CL:samtools view -H task3/data/tmp/vi0006.genome.markdup.downsampled.bam
# 4. Удаляем временную папку
shutil.rmtree(extracted_dir)
Вывод¶
Ура! Мы нашли полную команду запуска выравнивания, где было явно указано, что геном сборки hg38!
Решение¶
- Получение информации о покрытии и фильтрация - с помощью
bedtools genomecovполучаем.bedфайл с покрытием регионов, после чего отбираем только те, где покрытие >x10. - Определение координат APOE
- Определение сборки - с помощью
samtoolsв .bam файле находим информацию о сборке генома - Поиск координат APOE в соответствующей сборке генома - можно было просто зайти на NCBI, или его запарсить. Но! Есть прекрасная библиотека в питоне mygene, которая делает все сама. Вообще MyGene — это веб-тул, который построен поверх множества баз данных (NCBI, Ensembl, UniProt, ClinVar etc). Через REST API он отправляет запросы в эти БД по гену (или его символу, или NCBI/RefSeq идентификатору) и подтягивает нужную нам информацию (координаты для разных сборок, названия гена, варианты и связь с заболеваниями (в случае ClinVar и подобных БД), Gene Ontology и прочее). В общем, крутая штука! Горжусь братьями по делу!
- Определение сборки - с помощью
- Пересечение APOE и найденных регионов с покрытием >x10
- Создание
.bedфайла с координатами APOE - Поиск пересечений - с помощью
bedtools intersectвытаскиваем участки из нашего .bed файла, соответствующие APOE
- Создание
DATA_FOLDER_TASK4 = "task4/data"
task4_bampath = os.path.join(DATA_FOLDER_TASK4, "task4.bam")
Воспользуемся решением прошлого задания для определения сборки .bam-файла
bam_assembly_finder = BamAssemblyFinder(bam_file_path=task4_bampath)
_ = bam_assembly_finder.find_assembly()
Найдена сборка (по тегу @SQ в поле AS): NCBI37
Странная аннотация, но окей, похоже на GRCh37. Так что кажется наша сборка - hg19
P.S.: здесь тоже удивились такому странному id, потому что он не указывает, что это именно человека :)
Теперь наконец-то найдем долю гена APOE с покрытием >x10. Для этого воспользуемся GeneCoverageAnalyzer (подробнее в файле task4/analyze_gene_cov.py)
gene_cov_analyzer = GeneCoverageAnalyzer(bam_file_path=task4_bampath,
gene="APOE",
coverage_threshold=10,
genome_assembly="hg19",
species="human")
gene_cov_info = gene_cov_analyzer.compute(return_info=True)
Файл с координатами гена APOE уже существует по пути task4/data/APOE_human.bed Файл с отфильтрованными регионами по покрытию уже существует по пути task4/data/task4_x10.bed Файл пересечения уже существует по пути task4/data/task4_APOE_human.bed Длина гена APOE: 3639 п.о. Длина пересечения гена APOE и .bam-файла: 647 п.о. Доля гена APOE с покрытием больше 10x: 17.78%
gene_cov_info
{'intersected_length': 647,
'gene_length': 3639,
'coverage_fraction': 0.17779609782907393}
Я сразу оговорюсь - здесь доля покрытия может чуть варьироваться, потому что:
!cat {gene_cov_analyzer.gene_coord_file_path}
19 45409011 45412650
Вывод¶
Какое-то достаточно низкое покрытие.
Решение¶
В принципе схема аналогичная, как в задании №3. Открываем header .vcf, только уже с помощью bcftools, и ищем там информацию о сборке. Можно двумя путями:
- Ищем assembly напрямую (я видела это в информации о контигах)
- Лио ищем референс - навание файла, в котором может быть отражена сборка
DATA_FOLDER_TASK5 = "task5/data"
task5_vcfpath = os.path.join(DATA_FOLDER_TASK5, "task5.input.vcf.assembly.zip")
extracted_dir = os.path.join(DATA_FOLDER_TASK5, "tmp")
Эх и снова .zip. Делаем ту же схему, что и выше: распаковка $\rightarrow$ временная папка $\rightarrow$ делаем дела $\rightarrow$ удаляем временную папку
# Распаковка файла во временную папку
with zipfile.ZipFile(task5_vcfpath, 'r') as zip_ref:
zip_ref.extractall(extracted_dir)
print(f'Files: {os.listdir(extracted_dir)}')
Files: ['ki3640.imputed.vcf.gz']
# Указываем путь к файлу, находящимся во временной папке
task5_vcfpath = os.path.join(extracted_dir, "ki3640.imputed.vcf.gz")
# 3. Наконец-то найдем информацию о сборке
vcf_assembly_finder = VCFAssemblyFinder(vcf_file_path=task5_vcfpath)
vcf_assembly_finder.find_assembly_info()
Найдена сборка (по тегу assembly): hg19
И снова hg19!
Можем попробовать и по референсному файлу, но вообще это вариант, когда вариант с поиском assembly не сработает:
vcf_assembly_finder.find_reference()
Найдено референс (по тегу ##reference): file:///mnt/db/hg19/hg19.fa ##source="beagle.18May20.d20.jar" ##FILTER=<ID=LowQual,Description="Low quality"
'file:///mnt/db/hg19/hg19.fa\n##source="beagle.18May20.d20.jar"\n##FILTER=<ID=LowQual,Description="Low quality"'
Не забываем удалять временную папку:
shutil.rmtree(extracted_dir)
Решение¶
Для попарного сравнения генотипов будем использовать bcftools gtcheck. И здесь встает важный вопрос - запускаем с флагом -E 0 или без него?
- По умолчанию работаем с флагом
-E 40:- Существует ненулевая ($10^{-4}$) вероятность того, что несовпадение между генотипами вызвано ошибкой генотипирования.
- В таком случае Discordance может быть и дробным числом
- С флагом
-E 0:- Каждое несоответствие между генотипами считается жестким несовпадением.
- В таком случае Discordance может быть инетрпретирован как количество несовпадающих генотипов (ссылка), что число всегда целое.
Этому вопросу посвящено очень много обсуждений на гитхабе (#1538, #1556, #2194). И люди радовались, когда получали точное число несовпадающих генотипов (хотя отношение между Discordance для разных пар сохраняется). В общем, используем флаг -E 0*!
*Раньше -E был -e. В новой версии bcftools -e стал означать --exclude, так что нужно использовать тот или иной флаг в зависимости от своей версии bcftools.
**Кстати есть прикол, что bcftools gtcheck работает только с биаллельными сайтами, а мультиаллельные просто игнорирует (что обсуждалось в #1152 и сообщено в их коде). Можно предварительно "разбить" все мультиаллельные сайты на несколько строчек (bcftools norm -m-), как бы сделая их биаллельными. Я это прогоняла, и в наших образцах в принципе нет мультиаллельных, так что в нашем случае эта процедура избыточна :)
DATA_FOLDER_TASK6 = "task6/data"
task6_vcfpath = os.path.join(DATA_FOLDER_TASK6, "task6.vcf.gz")
comparator = VCFComparator(vcf_file_path=task6_vcfpath)
variant_comparison_df = comparator.compare(save_result=True)
Результат сравнения образцов из .vcf записан в task6/data/pairwise_comparisons.csv
variant_comparison_df
| Sample1 | Sample2 | Discordance | Average_neg_log_P_HWE | Number_of_sites_compared | Number_of_matching_genotypes | |
|---|---|---|---|---|---|---|
| 0 | HG00732 | HG00733 | 41781.0 | 0.030539 | 1066557 | 1024776 |
| 1 | NA20320 | HG00733 | 73481.0 | 0.014936 | 1066557 | 993076 |
| 2 | NA20320 | HG00732 | 74732.0 | 0.014677 | 1066557 | 991825 |
Что эти колонки вообще обозначают?¶
- Sample1
Первый образец (его идентификатор) из сравниваемой пары
- Sample2
Второй образец (его идентификатор) из той же пары, с которым сравнивается первый образец
- Discordance
Метрика, отражающая, насколько сильно отличаются генотипы в паре образцов. В нашем случае, а именно при использовании флага -E 0, Discordance - это количество несовпадающих генотипов.
- Average_neg_log_P_HWE
Вероятность отклонения (а точнее минус логарифм от нее) генотипов от равновесия Харди-Вайнберга (HWE). Согласно HWE, в идеальной популяции, где спаривание происходит случайным образом и отсутствуют такие факторы как мутация, миграция, отбор и генетический дрейф, частоты генотипов остаются стабильными от поколения к поколению.
Хотя эти предположения редко бывают верны в естественном мире, они позволяют нам вычислить ожидаемую частоту аллеля. Значительные различия между наблюдаемыми и ожидаемыми частотами указывают на то, что «что-то» (т. е. одно или несколько из вышеперечисленного) происходит .
В контексте результата bcftools, для каждой позиции (локуса) в геноме у нас есть частота каждого аллеля и частота генотипов в этой позиции. На основе этих частот можно рассчитать ожидаемые частоты генотипов при условии HWE. Затем наблюдаемые частоты сравниваются с ожидаемыми, и на основе этого сравнения для каждой позиции оценивается P(HWE) - вероятность того, что частоты генотипов, наблюдаемые для этого сайта, соответствуют равновесию Харди-Вайнберга.
В bcftools, когда выполняется проверка на соответствие равновесию Харди-Вайнберга, используется $\chi^2$-тест:
$$\chi^2 = \sum{\frac{(O_i - E_i)^2}{E_i}}$$
где $O_i$ - наблюдаемые частоты, а $E_i$ - ожидаемые частоты
- Number_of_sites_compared
Общее количество генетических позиций (сайтов), по которым было выполнено сравнение между двумя образцами
- Number_of_matching_genotypes
Количество позиций, по которым генотипы образцов совпадают - сколько генетических позиций идентичны между двумя сравниваемыми образцами. В случае нашего флага -E 0:
$$\text{[\text{Number of sites compared}]} = \text{[Number of matching genotypes]} + \text{[Discordance]}$$
На этом моменте я искренне перестала понимать, почему люди с обсуждений на гите так счастливы использовать флаг -E 0, если количество несовпадающих генотипов можно и самому рассчитать... а с тем же дефолтным -E 40 можно хотя бы регулировать вероятность ошибки генотипирования
variant_comparison_df['Discordance'] + \
variant_comparison_df['Number_of_matching_genotypes'] \
== variant_comparison_df['Number_of_sites_compared']
0 True 1 True 2 True dtype: bool
Визуализация результата¶
1. Относительное число несоответствий¶
Нарисуем тепловую карту (heatmap) по Discordance. Более "яркие" значения соответствуют большим значениям Discordance - то есть такая пара образцов имеет большое количество несовпадающих генотипов, и наоборот. Самый маленький Discordance у пары образцов сам-с-собой (sample vs sample) и будет равным 0
# Получаем уникальные пары оразцов
samples = pd.unique(variant_comparison_df[['Sample1', 'Sample2']].values.ravel('K'))
# Создаем пустую симметричную матрицу с нулями
matrix = pd.DataFrame(np.zeros((len(samples), len(samples))), index=samples, columns=samples)
# Заполняем матрицу значениями Discordance
for _, row in variant_comparison_df.iterrows():
matrix.at[row['Sample1'], row['Sample2']] = row['Discordance']
matrix.at[row['Sample2'], row['Sample1']] = row['Discordance']
# И еще нормализуем на общее число сравниваемых сайтов - оно одинаково между образцами
# Идея такая, что хочется посмотреть на % несовпадающих сайтов, а не их абсолютное число
assert variant_comparison_df['Number_of_sites_compared'].nunique() == 1, \
"'Number_of_sites_compared' различается между образцами"
n_sites = variant_comparison_df['Number_of_sites_compared'].iloc[0].item()
matrix = 100 * matrix / n_sites
ax = sns.heatmap(matrix,
cmap = 'Purples', linewidths=1, linecolor="grey",
annot=True, fmt=".2f")
ax.set_xticklabels(matrix.columns, rotation=45)
ax.set_yticklabels(matrix.index, rotation=0)
ax.set_title("Относительное число несовпадений (%)")
plt.show()
Самое большое число несоответствий у пары NA20320-HG00732, и примерно такое же число несоответствий - у пары NA20320-HG00733. Конечно наша выборка только из трех образцов (и $C_3^2=3$ пар образцов), но что-то мне подсказывает, что NA20320 сильно отбивается от HG00732 и HG00733. Точнее не что-то, а 1) названия образцов, 2) значения Discordance
2. Соблюдение HWE¶
Не будем отклоняться от прошлого графика - также нарисуем тепловую карту средних значений вероятности отклонения от HWE
# Создаем пустую симметричную матрицу для Average_neg_log_P_HWE
hwe_matrix = pd.DataFrame(np.zeros((len(samples), len(samples))), index=samples, columns=samples)
# Заполняем матрицу значениями Average_neg_log_P_HWE
for _, row in variant_comparison_df.iterrows():
hwe_matrix.at[row['Sample1'], row['Sample2']] = row['Average_neg_log_P_HWE']
hwe_matrix.at[row['Sample2'], row['Sample1']] = row['Average_neg_log_P_HWE']
# Визуализируем тепловую карту
ax = sns.heatmap(hwe_matrix,
cmap='Blues', linewidths=1, linecolor="grey",
annot=True, fmt=".2f")
ax.set_xticklabels(hwe_matrix.columns, rotation=45)
ax.set_yticklabels(hwe_matrix.index, rotation=0)
ax.set_title(r"AVG $\left[ -\log{P(\text{HWE})} \right]$", fontsize=14)
plt.show()
Чтобы интерпретировать картинку, сделаем небольшое вступление (на всякий случай, я сама порой путаюсь):
$P(\text{HWE})$ - вероятность отклонения от HWE для пары образцов.
- Бóльшая вероятность указывает на соблюдение HWE у двух образцов по этой позиции
- Мéньшая вероятность указывает на отклонение от HWE у двух образцов по этой позиции
$-\log{P(\text{HWE})}$ - то же самое, что и выше, только с мат. наворочками над исходной вероятностью:
- Низкие значения соответствуют более высокому p-value, что указывает на то, что нет значимых отклонений от HWE у двух образцов по этой позиции.
- Высокие значения соответствуют более маленькому p-value, что указывает на то, что с какой-то вероятностью можно говорить о значимых отклонениях от HWE в этой паре образцов по данной позиции.
На нашей тепловой карте мы видим только низкие значения -$\log{P(\text{HWE})}$, что указывает на отсутствие значительных отклонений от равновесия Харди-Вайнберга среди всех пар образцов. А если у нас нет значительных нарушений HWE, то это может указывать на отсутствие значительных генетических различий или мутаций между образцами в сравниваемых позициях.
Выводы¶
Мы посмотрели на красивые картинки и сделали два отдельных вывода:
- Для пар с образцом NA20320 относительное число несовпадений составляет примерно 7%, тогда как для последней пары — около 4%. Эти значения могут показаться значительными на фоне числа сравниваемых позиций, которое составляет порядка миллиона
- Все пары образцов показывают низкие значения $-\log{P(\text{HWE})}$, что указывает на незначительные отклонения от равновесия HW в рассматриваемых позициях
Таким образом, несмотря на высокое абсолютное и относительное количество несовпадающих позиций между образцами, низкие значения $-\log{P(\text{HWE})}$ свидетельствуют о том, что эти несовпадения не являются статистически значимыми. То есть различия в генотипах не указывают на значимые генетические изменения или мутации, а скорее связаны с естественным генетическим разнообразием.
В общем, данные свидетельствуют о наличии генетического разнообразия между образцами, но не о значительных мутациях, что может указывать на близкое или умеренное родство между ними.
Решение¶
DATA_FOLDER_TASK7 = "task7/data"
task7_path = os.path.join(DATA_FOLDER_TASK7, "task7.zip")
extracted_dir = os.path.join(DATA_FOLDER_TASK7, "tmp")
Мы снова имеем дело с .zip-архивом, поэтому следуем той же логике, что и ранее: распакуем его во временную папку, сделаем дело, и после временную папку подчистим. Я написала кастом для сохранения таблицы частотности аллелей, поэтому результат можем сохранить. Но промежуточное точно не нужно
with zipfile.ZipFile(task7_path, 'r') as zip_ref:
zip_ref.extractall(extracted_dir)
print(f'Files: {os.listdir(extracted_dir)}')
Files: ['data']
О боже, папка в папке.. Ну ладно, сделаем ход конем и вручную переместим файлы в tmp/
!mv {os.path.join(extracted_dir, "data/*")} {extracted_dir}
!rm -r {os.path.join(extracted_dir, "data/")}
Супер! Теперь, наконец, определяем пути к файлам. Только перед этим посмотрим на них:
os.listdir(extracted_dir)
['ALL.chr15.subset.hg38.annotated.clean.vcf.gz', 'tubes_pop.tsv']
population_data_path = os.path.join(extracted_dir, 'tubes_pop.tsv')
frequency_data_path = os.path.join(extracted_dir, 'ALL.chr15.subset.hg38.annotated.clean.vcf.gz')
af_savepath = os.path.join(DATA_FOLDER_TASK7, 'AF_populations.csv')
И приступим к самому главному - подсчитаем частоты! (под капотом в основном bcftools)
genvar_analyzer = GeneticVariantAnalyzer(population_data_path=population_data_path,
frequency_data_path=frequency_data_path)
Количество образцов в анализируемых популяциях: pop AFR 671 EUR 521 SAS 492 Name: count, dtype: int64
AF_frequencies = genvar_analyzer.calculate_af(save_path=af_savepath)
Подсчет генетических частот...: 100%|██████████| 3/3 [05:35<00:00, 111.67s/it]
# посмотрим, как выглядит итоговая табличка
AF_frequencies.head(5)
| CHROM | POS | REF | ALT | AF | Population | |
|---|---|---|---|---|---|---|
| 0 | 15.0 | 17009562.0 | T | C | 0.000000 | EUR |
| 1 | 15.0 | 17009590.0 | G | T | 0.445298 | EUR |
| 2 | 15.0 | 17009613.0 | A | C | 0.000000 | EUR |
| 3 | 15.0 | 17009614.0 | T | C | 0.000960 | EUR |
| 4 | 15.0 | 17009622.0 | G | A | 0.000000 | EUR |
Перед тем, как мы предпримем попытку что-то нарисовать и интерпретировать, посмотрим, что у нас вообще за популяции. И здесь интересно, зачем образцы из SGDP отделили... возможно, чтобы контролировать батч-эффект, но это не точно
genvar_analyzer.population_data[['pop', 'pop_name']].value_counts()
pop pop_name AFR African Ancestry 663 EUR European Ancestry 514 SAS South Asian Ancestry 486 AFR African Ancestry,Africa (SGDP) 8 EUR European Ancestry,West Eurasia (SGDP) 7 SAS South Asia (SGDP),South Asian Ancestry 6 Name: count, dtype: int64
График рассеяния (scatter plot)¶
genvar_analyzer.plot_scatter(df_frequencies=AF_frequencies, populations=genvar_analyzer.populations)
красота 👍
Что за точки и как на это смотреть?¶
Думаю не будет вредно разобраться, что означает каждое положение точки на графике. Из формулировки задания освежим, что точка - это генетический вариант с соответствующими частотами (частоты варианта в популяции, отраженной на осях OX и OY). Цветом обозначена позиция (чем темнее - тем ближе к началу 15-й хромосомы, чем светлее - тем ближе к концу 15-й хромосомы). Но так как вариантов очень много, последние точки (т.е. расположенные ближе к концу 15-й хромосомы, и, соответственно более желтые) легли на графике на все остальные. Тем не менее по краям мы можем увидеть точки, соответствующие середине или началу хромосомы :)
Итак, к графику (этот текст я писала, когда я неправильно доставала генетические частоты, и график был более красивый, но сейчас тоже +- актуально):
Точка лежит на линии регрессии или близка к ней, т.е. $y=x$. Если точка находится на этой линии, то в обеих популяциях ~одинаковая распространненность данного генетического варианта. То есть вариант:
- Возможно имеет стабильное распределение вне зависимости от популяции
- Возможно нейтральная мутация, не подвергающаяся селекционному давлению ни в одной из популяций внутри пары.
Точка сильно выше линии $y=x$. Тогда генетический вариант более распространен во второй популяции (OY) по сравнению с первой (OX). То есть вариант:
- Возможно более благоприятен в среде второй популяции, где он дает преимущество (в то время как во второй популяции данный вариант или вреден, или нейтрален $\rightarrow$ меньше частота варианта в этой популяции)
- Скорее всего вариант имеет важное биологическое значение и подвергается разному давлению выбора в этой паре популяций
Точка сильно ниже линии $y=x$. Тогда генетический вариант более распространен в первой популяции (OX) по сравнению со второй (OY). В принципе схема та же, что и в п.2
Выводы по графику рассеяния¶
- Если бы мы построили график частот аллелей в одной популяции против частот тех же аллелей в этой же популяции, мы бы получили прямую линию ($y=x$). Соответственно, графики, которые больше прижаты к линии $y=x$, указывают на более высокую генетическую схожесть между сравниваемыми популяциями
- График рассеяния для пары популяций SAS-EUR наиболее близок к линии $y=x$. Так что, вероятно, популяции SAS и EUR являются наиболее генетически схожими из всех (трех) исследуемых пар популяций
- Графики для пар популяций AFR-EUR и AFR-SAS занимают значительную часть графика, почти заполняя квадрат [0, 0] - [1, 1]. То есть между этими парами популяций существует значительное различие в частотах генетических вариантов. При этом в обеих парах наблюдается как высокая частота аллелей в популяции AFR и низкая в другой (SAS/EUR), так и наоборот. Скорее всего популяция AFR значительно отличается от популяций EUR и SAS, возможно из-за сильных различий в среде обитания
3D график рассеяния¶
Мне просто стало интересно...
Если не запускаете ноутбук (который, возможно, и не запустится из-за зависимостей, которые на другой машине не могу проверить + некоторые файлы я переименовывала), то можете открыть и поиграться с интерактивными картинками, которые я сохранила в формате .html. Они находятся по этой ссылке
Дальше есть два варианта действий:
- Скачать любой файл и открыть локально. Но не факт, что откроется с телефона
- Добавить приставку https://html-preview.github.io/?url= к ссылке на .html` Например: https://html-preview.github.io/?url=https://github.com/ksumarshmallow/Practical-BF/blob/main/HW1/task7/imgs/AF_3D_17009562_25506717.html
# Собираем частотные таблицы по каждой популяции
dfs = genvar_analyzer._split_and_filter_dfs(df_frequencies=AF_frequencies, populations=genvar_analyzer.populations)
# Названия популяций
pops = [k for k in dfs.keys()]
# AF для каждой популяции
AFs = np.array([df.AF.values for df in dfs.values()])
# Позиции
positions = dfs[pops[0]].POS.values.astype(int)
AFs.shape
(3, 2294447)
Все позиции мы не сможем нарисовать (все сломается и будет некрасиво), поэтому нарисуем несколько графиков!
# Создаем папку для сохранения .html, если ее еще нет
AF_3d_folder = os.path.join('task7', 'imgs')
if not os.path.exists(AF_3d_folder):
os.makedirs(AF_3d_folder)
Наконец-то рисуем интерактив! Будем идти с шагом 100.000 по генетическим вариантам
step = 100_000
kwargs = {"x_title": pops[0], "y_title": pops[1], "z_title": pops[2]}
for i in tqdm(range(0, len(positions), step), desc="Plot 3D AFs...", colour="CYAN"):
x = AFs[0][i:i+step]
y = AFs[1][i:i+step]
z = AFs[2][i:i+step]
pos = positions[i:i+step]
plot_name = f'AF_3D_{pos[0]}_{pos[-1]}'
# рисуем график
fig = plot_3d_scatter(x=x, y=y, z=z, c=pos, plot_name= plot_name, **kwargs)
# сохраняем в папочку
fig.write_html(f"{AF_3d_folder}/{plot_name}.html")
Plot 3D AFs...: 100%|██████████| 23/23 [00:03<00:00, 6.60it/s]
Посмотрим на последнюю картинку, делая сначала ноутук интерактивным:
- В другой среде возможно нужно будет поставить
%matplotlib widget - Для просмотра с гитхаба инструкция выше
%matplotlib notebook
fig
Особенности на 3D графике рассеяния¶
- Лучше на него смотреть, ориентивовав "на себя" главную диагональ
- Можно регулировать график - наверху для этого есть виджеты (Zoom, Pan=перемещать сам кубик, Orbital rotation, Turntable rotation)
- Если открыли в интерактиве - можно навести курсор на точку и посмотреть, какие частоты соответствуют этому варианту у всех трех популяций
- Отклонения от диагонали - варианты, которые отличаются как минимум между двумя популяциями. Они нам наиболее интересны!
Выбросы (outliers)¶
Меня очень сильно понесло (ночное вдохновение) и сейчас я хочу посмотреть на "выбросы" - те генетические варианты, которые отбиваются от остальной кучи точек на графике рассеяния и 3D графике рассеяния.
В целом выбросы интересны потому, что они являются вариантами, где частота аллеля сильно отличается между популяциями. То есть это что-то индивидуальное у одной популяции и может быть следствием особых условий среды в этой популяции. Вообще мне легче будет найти выбросы из 3D граффике рассеяния. Поэтому формулирую следующее:
Цель данного раздела - найти генетические варианты, частота которых различается хотя бы между двумя популяциями (из EUR, AFR, SAS).
В контексте выбросов по генетическим вариантам самый популярный пример - вариант, связанный с геном HBB, который кодирует β-цепь гемоглобина. Мутации в последнем:
- Приводят к различным заболеваниям крови, включая серповидноклеточную анемию;
- Оказывает сильное влияние на восприимчивость к малярии.
То есть, существует связь между серповидно-клеточной анемей и защитой от малярии. Этот вариант особенно распространен среди жителей Западной и Центральной Африки, некоторых районов Индии и Средиземноморья, где малярия была исторически распространена. Люди с одной копией мутантного аллеля (гетерозиготы) имеют частичную защиту от малярии, что делает этот вариант частым в популяциях, подвергшихся сильному давлению естественного отбора.
Так что это типичный пример того, как различия в окружающей среде могут приводить к резкому расхождению частот аллелей между популяциями. И я хочу посмотреть, что же происходит у нас!
Выбросы заметны на построенном 3D scatter-plot'е, но их нам нужно как-то вычислить. Я пойду наивным образом:
- Извлечем частоты для каждого генетического варианта в каждой популяции
- Вычислим дисперсию генетических частот среди популяций (да, по трем, но а что делать)
- Генетические варианты с самой высокой дисперсией - те самые варианты, частоты которых сильно различаются между тремя популяциями
# Первый шаг - сделан выше
# Второй шаг - считаем стандартное отклонение
deviations = AFs.std(axis=0)
print(f'Максимальная дисперсия: {deviations.max():.4f}')
# Третий шаг:
# 1. Cчитаем порог по отклонению
dev_threshold = deviations.max() - deviations.std() # np.quantile(deviations, 0.995).item()
print(f'Порог по дисперсии: {dev_threshold:.4f}')
# 2. Определяем позиции, где отклонение > порога
deviations_filtered_idx = np.where(deviations > dev_threshold)[0]
positions_filtered = positions[deviations_filtered_idx]
print(f'Количество генетических вариантов после отсечки: {positions_filtered.shape[0]}')
# 3. Находим в них соответствующие частоты
AFs_filtered = AFs[:, deviations_filtered_idx]
Максимальная дисперсия: 0.4175 Порог по дисперсии: 0.3847 Количество генетических вариантов после отсечки: 11
Выведем наши самые крутые позиции!
outliers_variants = pd.DataFrame(AFs_filtered, index = pops, columns = [f'AF_{i+1}' for i in range(len(positions_filtered))]).T
outliers_variants['chr15_pos'] = positions_filtered
outliers_variants
| EUR | SAS | AFR | chr15_pos | |
|---|---|---|---|---|
| AF_1 | 0.050864 | 0.019309 | 0.850969 | 25275198 |
| AF_2 | 0.048944 | 0.019309 | 0.851714 | 25275760 |
| AF_3 | 0.008637 | 0.001016 | 0.890462 | 29135197 |
| AF_4 | 0.972169 | 0.938008 | 0.119225 | 33966633 |
| AF_5 | 0.972169 | 0.938008 | 0.131893 | 33967969 |
| AF_6 | 0.976008 | 0.940041 | 0.137854 | 33973629 |
| AF_7 | 0.976967 | 0.940041 | 0.142325 | 33977774 |
| AF_8 | 0.976967 | 0.958333 | 0.151267 | 33979447 |
| AF_9 | 0.935701 | 0.966463 | 0.102832 | 45095352 |
| AF_10 | 0.926104 | 0.960366 | 0.093890 | 45105368 |
| AF_11 | 0.077735 | 0.115854 | 0.913562 | 60563356 |
Ладно, мой способ нахождения был не самым лучшим, потому чтос самая высокая дисперсия будет там, где хотя бы в паре популяций у одной частота варианта будет >0.9, а у второй <0.1. Но все равно мы нашли, по-моему, что-то интересное! То есть именно те варианты, которые часто встречаются в одной популяции, но редко встречаются в других!. Либо наоборот - редко встречаются в одной популяции, но часто бывают в других. И, кстати, снова от других популяций отличается только AFR
Теперь мы бежим в NCBI Genome Data Viewer и смотрим на наши позиции - возможно, для них уже что-то известно. Для всех нашлись rsid. И табличка по AF среди популяций в любых случаях совпадает с нашей!
Я все обобщила тут:
| Название гена | Позиция | rsID и ссылка | Тип (SNV/индел) | Публикации | ClinVar |
|---|---|---|---|---|---|
| SNHG14 | chr15:25275198 | rs1402074614, rs1442810789, rs2073056425, rs577002978 | Индел, SNV | Нет публикаций | Нет |
| chr15:25275760 | rs2714786 | SNV | Нет публикаций | Нет | |
| chr15:29135197 | rs10152250 | SNV | Нет публикаций | Нет | |
| chr15:33966633 | rs600637 | SNV | Нет публикаций | Нет | |
| AVEN, CHRM5 2KB Upstream Variant | chr15:33967969 | rs676281 | SNV | Нет публикаций | Нет |
| chr15:33973629 | rs1175001719, rs476814 | Индел, SNV | Нет публикаций | Нет | |
| chr15:33977774 | rs480616 | SNV | Нет публикаций | Нет | |
| chr15:33979447 | rs592585 | SNV | Нет публикаций | Нет | |
| DUOX2 | chr15:45095352 | rs199138 | SNV |
1. A Population-Based Genomic Study of Inherited Metabolic Diseases Detected Through Newborn Screening 2. Effective selection of informative SNPs and classification on the HapMap genotype data |
ClinVar (rs199138) |
| chr15:45105368 | rs269863 | SNV | Нет публикаций | ClinVar (rs269863) | |
| RORA, RORA-AS1 | chr15:60563356 | rs974828 | SNV | Нет публикаций | Нет |
Очищаем пространство¶
# Удаляем временную папку
shutil.rmtree(extracted_dir)
Вывод¶
- Мы научились вычислять частоты генетических вариантов для различных популяций
- Мы визуализировали частоты на графике рассеяния и 3D графике рассеяния (и поигрались с последним)
- Из визуализации частот мы поняли, что популяция AFR отличается по частоте генетических вариантов от популяций EUR и SAS
- Мы нашли наиболее варьирующиеся позиции, которые я назвала выбросами - те генетические варианты, которые стоят далеко от главной диагонали на 3D графике рассеяния.
- Каждую позицию мы просмотрели в NCBI DATA Viewer, определив, в каком гене находится данный генетический вариант, и найдя соответствующий rsid
- Еще раз убедились, что популяция AFR более всего отличается от других двух популяций
Итог¶
Я исчерпала все свое вдохновение и силы программировать, но я в восторге! 🌟